/*
 * This class represents a number of combinations as integers, longs or some other
 * objects that can be converted from integers. In this case, there may be initialized a set
 * of known factorials. Otherwise, they must be calculated initially. This feature should be
 * achieved with the aid of inheritance or some macro design
 */
#ifndef _COMBINATION_H
#define _COMBINATION_H

template<typename T = int> 
class combination{
      public:
             combination(T k, T l):n(k),r(l){
                    facVec[0] = T(1);facVec[1] = T(1);facVec[2] = T(2);facVec[3] = T(6);facVec[4] = T(24);
                    facVec[5] = T(120);facVec[6] = T(720);facVec[7] = T(5040);facVec[8] = T(40320);facVec[9] = T(362880);
                    facVec[10] = T(3628800);facVec[11] = T(39916800);facVec[12] = T(479001600);       
                    result = evaluate(k,l);         
             }
             T get_result();  //Return the value of combination
      private:
              T n;  //Total number
              T r;  //Number to combine
              static const int knownnumb = 13;  //Number of known factorials
              T facVec[knownnumb];
              T result;  //To store the result for this combination class
              T evaluate(T a, T b);
};

//Interface function for combination class
template<typename T>
T combination<T>::get_result(){
    return result;
}
//Core funtion beneath combination class
template<typename T>
T combination<T>::evaluate(T a, T b){
    if(a<b) return T(0);
    if(a<T(0) || b<T(0)) return T(0);
    if(b==T(0) || a==b) return T(1);
    if(b==T(1) || b == a-T(1)) return a;
    if(T(2)*b>a) b = a-b;
    if(b<T(knownnumb)){
        T res = T(1);
        int k;
        for(k=0;T(k)<b;k++){
                res = res *(a-T(k));
        }
        return res / T(facVec[k]);
    }
    //Recursive calculations
    T p, q, temp;
    p = b / T(2);
    q = b - p;
    temp = a - p -q + T(1);
    T rr = T(0);
    for(int i=1; T(i)<temp;i++){
            rr = rr + evaluate(temp+p-T(i)-T(1),p-T(1))*evaluate(T(i)+q-T(1),q);
    }
    return rr;
}
#endif
